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Nambu-Goldstone modes associated with (topological) defects such as vortices and domain walls in (su- 
per)fluids are known to possess quadratic/non-integer dispersion relations in finite/infinite-size systems. Here, 
we report interpolating formulas connecting the dispersion relations in finite- and infinite-size systems for Kelvin 
modes along a quantum vortex and ripplons on a domain wall in superfluids. Our method can provide not only 
the dispersion relations but also the explicit forms of quasiparticle wavefunctions ( u , v). We find a complete 
agreement between the analytical formulas and numerical simulations. All these formulas are derived in a fully 
analytical way, and hence not empirical ones. We also discuss common structures in the derivation of these 
formulas and speculate on the general procedure. 
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I. INTRODUCTION 

In the latter part of the 19th century. Lord Kelvin left many 
influential works in classical fluid mechanics, and not a few of 
them form the foundation of this research field today. Among 
his works are those on the propagation of linear waves in the 
vicinity of local inhomogeneous structure of fluids, such as 
ripple modes (capillary waves) along an interface between two 
fluids 1 , which arise as a by-product of the study of Kelvin- 
Helmholtz instability (e.g.. Ref. 2), and helical motions of 
vortices, which are now called Kelvin modes 3 . These modes 
are notable for the point that they have non-integer dispersion 
relations: while the ripple modes have a fractional dispersion 
relation e oc /c 3/2 (Ref. 4), the Kelvin modes have a logarith¬ 
mic one e oc -kr log k. 

In modern physics, these linear waves are also known to 
emerge in various examples of quantum fluids. The Kelvin 
modes, or Kelvons if observed as quantized quasiparticles, 
exist in quantized vortices in superfluids 5-8 , Bose-Einstein 
condensates (BECs) of ultracold atomic gases 9-11 , or neu¬ 
tron superfluids in neutron stars, having the same disper¬ 
sion relation with classical fluids in the infinite-volume limit. 
Kelvin modes are considered to play an important role known 
as the Kelvin-mode cascade in turbulences, including quan¬ 
tum turbulence 1213 . Thus, understanding Kelvin modes bet¬ 
ter is an important step toward complete characterization of 
turbulences, which remains an unsolved problem since the 
first observation by da Vinci. The ripple modes, or rip¬ 
plons if identified as quasiparticles, emerge on a domain 
wall 14-16 (DW) of a mixture of two kinds of BECs and also 
possess the same dispersion relation with classical fluids in 
infinite-size systems 17,18 , and the analogous phenomena of the 
Kelvin-Helmholtz and Rayleigh-Taylor instabilities were also 
found 19-21 . There are also related issues 22-24 . 

Recently, a new insight has been brought to these gapless 
modes, stimulated by a renewed understanding on Nambu- 
Goldstone modes (NGMs) in non-relativistic systems 25-30 . 


Both Kelvin modes 31,32 and ripple modes 18,32 have quadratic 
dispersion relations, e ~ (log R)k 2 and e ~ sfZk 2 , in finite-size 
systems with R and L denoting system lengths perpendicular 
to a vortex and DW, respectively. These facts are consistent 
with the general argument that an NGM with quadratic disper¬ 
sion corresponds to two broken symmetries 26-29 . In the limit 
R, L —> oo, however, we encounter a difficulty of the diver¬ 
gent coefficient and the correct dispersion laws change to the 
non-integer ones mentioned above. How these qualitatively 
different integer and non-integer laws are continuously inter¬ 
polated is yet to be clarified. The finite-size correction will be 
also crucial for quantum turbulences with a large number of 
vortices, since the mean intervortex distance gives the effec¬ 
tive system size for each vortex. 

In this paper, we report analytical formulas interpolating 
the integer and non-integer dispersions in finite- and infinite- 
size systems for Kelvin and ripple modes, and find a complete 
agreement with numerical simulations. We also summarize 
common practical procedures in derivation of these two ex¬ 
amples, which could become a guiding principle to derive in¬ 
terpolating formulas for NGMs around other topological de¬ 
fects. 

The organization of this paper is as follows. In Sec. II, we 
summarize our main analytical formulas and their numerical 
verifications for Kelvin modes and ripplons. We also sum¬ 
marize common aspects of mathematical derivations given in 
subsequent sections. In Secs. Ill and IV, we provide full 
analytical derivations of main results for Kelvin modes and 
ripplons, respectively. Section V is devoted to a summary. 
Appendices A and B provide a few technical calculations for 
DWs in two-component BECs. 
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II. MAIN RESULT AND NUMERICAL EVIDENCE 
A. Kelvin modes 


First we report the interpolating dispersion formula for 
Kelvin modes propagating along a quantized vortex. The de¬ 
tailed derivations are given in Sec. III. We consider an in¬ 
finitely long cylinder with radius R. The Gross-Pitaevskii 
(GP) energy functional for a single-component BEC with 
chemical potential term is given by 

H -, N = J d i r ^ +gW *-MA. (2.1) 

Without loss of generality we set 2m = j = g = 1 by 
rescaling of variables. The GP equation is then given by 
id,if/ = -V 2 / - 2iA + 2\i//\ 2 ifs. The boundary condition (BC) at 
r = R does not affect the main results shown below. For exam¬ 
ple, it can be either Dirichlet or Neumann. We are interested 
in a stationary single vortex solution. Setting <// = /(r)e lfl , the 
function / satisfies -f" - 7 - + ^ - 2/(1 - f 2 ) = 0. Hence¬ 
forth, we write the vortex solution in the infinite-size system 
(R = 00 ) as f m (r). The asymptotic form for large r is given by 
foo(r) = 1 - -jk + ()(r 4 ). The Bogoliubov equation 33-36 de¬ 
scribing quasiparticle excitations is obtained by substituting 
/ = / + z<e~ lc ' + v*e lf *' into the GP equation and lineariz¬ 
ing it with respect to ( u , v). Then, our main result for Kelvin 
modes is summarized as follows. The dispersion relation e k 
and the quasiparticle wavefunctions for Kelvin modes, which 
we write (n, v) = (u k (r), v k (r)e~ 2 ' e )e lkzZ , are given by 

( 2 . 2 ) 
(2-3) 

K 0 (k)+K 2 (k) 
h(k)+I 2 (k) ’ 

(2.4) 


e k =k 2 (- log | + ri - y-xikR)) , 

WL/ftW-l + f+/») 

M r V \-Fk(r) + \-^+fL(r))' 

F k {r) :=k[Ki(kr) + x(kR)h(kr)], x( k ) := 


where k - \k z \, I n ,K n are the modified Bessel function of the 
first and second kind, y = 0.577 ... is the Euler-Mascheroni 
constant, and 77 is a constant defined by 


X oo 

dr [r /» 2 - 2/<x>(r)//(r) log r] - 0.227. (2.5) 

Since x(k) has the expansion 

fir + (~y - i -logl) + 0 (k 2 ) (k « 1 ) 


*(*) 


- ) P 


-2k 


[1 + 


i+o(k- 2 )] 


(k »i), 


(2.6) 


the dispersion formula [Eq. (2.2)] includes the following two 
important limiting cases: 

^ f-#+ k 2 (logR + | +rj) (kR « 1) (2.7a) 
6k ~ I k 2 (- log | + 77 - y) (R-> 00 ). (2.7b) 

The expression (2.7a) revisits the result of Refs. 31 and 32, 
except for the correction term | + 77 for the //-coefficient. 



FIG. 1. (Color online) The dispersion e k of Kelvin modes for R = 
10 and R = 100 under the Neumann BC. The analytical formula 
[Eq. (2.2)] and numerical solutions agree well in the low k region. 
Equations (2.2) with R = 100 and Eq. (2.7b) for R —> 00 are almost 
the same and the two lines for them in the figure overlap each other. 


The expression (2.7b) describes the non-integer dispersion 
e ~ -k 2 log k in the infinite volume 5-8 . The correction terms 
including 77 improve the fitting with numerical results. This 
constant is slightly different from the previously-known value 
2 (Ref. 6 ); this difference arises from the use of explicit quasi¬ 
particle wavefunctions Eq. (2.3). The equivalent expression 
for this 77 was also reported in Ref. 37. The formula (2.2) well 
explains numerical data not only for the above-mentioned lim¬ 
iting cases but also for the intermediate regions. See Fig. 1. 

The quasiparticle eigenstate [Eq. (2.3)] with R — 00 in¬ 
cludes Pitaevskii’s result 5 in two ways; First, setting k — 0, it 
reduces to (u 0 (r),v 0 (r)) = (// + ~ 7 O’ whic h has the 

physical meaning of the zero-mode solution originated from 
translational symmetry breaking 32 . (See also Sec. Ill A of this 
paper.) Second, if we focus on the asymptotic region r » 1, 
we have (u k (r),v k (r)) oc K\(kr), which was used to derive 
e - k 2 log k in Ref. 5. While F k (r) has a power series with re¬ 

spect to k if R < 00 , it becomes invalid for R - 00 , since K\ ( kr ) 
has a logarithmic term. This means that the naive perturbative 
expansion does not work when R = 00 . Equation (2.3) well 
explains the numerical solutions for quasiparticle excitations. 
See Fig. 2. 

While the numerical results shown in Figs. 1 and 2 are 
those under the Neumann BC [i.e., f'(R) - u'(R) - v’(R) — 
0 ], our analytical results are also well applicable for the sys¬ 
tems obeying the Dirichlet BC [i.e., f{R) = u{R ) = v(R) = 0]. 
Analytical formulas without any modification can show a 
modestly good agreement with numerical results even for the 
Dirichlet BC. As we will see below, however, if we introduce 
an effective system radius R-/3 with a numerical fitting param¬ 
eter [3 - 0.946, we obtain a more refined agreement between 
the numerical results and the analytical formulas. 

Figure 3 shows the A’-dcpendence of the energy of zero- 
wavenumber solution f (l . For the Neumann BC, it is well fit¬ 
ted by the formula eo = —jp, consistent with Eq. (2.7a) and 
Ref. 31. For the Dirichlet BC, if we fit the numerical result 
by the ansatz — ^ 2 . we find 78 - 0.946. The physical mean¬ 
ing of this p is obvious; since the Dirichlet BC suppresses 
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FIG. 2. (Color online) The zero- and finite-wavenumber (k = 0.1) 
solutions of Kelvin modes under the Neumann BC. Here we set R = 
10 . 


FIG. 4. (Color online) The dispersion relation ek of Kelvin modes 
for R = 10 and R = 100 under the Dirichlet BC. Here, when we plot 
the analytical formulas Eqs. (2.2) and (2.7a), we use the modified 
system radius R —/? instead of the bare R. 


the wavefunctions near the boundary, the effective radius gets 
shorter than that of the Neumann BC by a length about the 
healing length. Figure 4 shows the comparison of dispersion 
relations between the numerical results and the analytical for¬ 
mulas with R being replaced by R-/3. The fitting is improved 
drastically by using R-/3 instead of the bare R. Figure 5 shows 
the quasiparticle wavefunctions, showing a good agreement 
with the analytical formulas except near the boundary. 


B. Ripplons 


Next, we report the dispersion relation of ripplons on a DW 
in two-component BECs. The details of the derivation are 
given in Sec. IV. The energy functional is given by 


H = 



+ X giM 2 m 2 

i,j= 1,2 


(2.8) 


Here we assume gn,gzi > 0 and g n = g 2 i > yjgugn, in 
which case the ground state is given by the state such that </q 



FIG. 5. (Color online) The zero- and finite-wavenumber (k = 0.1) 
solutions of Kelvin modes under the Dirichlet BC with the system 
radius R = 10. The analytical formula [Eq. (2.3)] is used with re¬ 
placing R by R - /I, where /? = 0.946 is obtained from the fitting in 
Fig. 3. Numerical solutions almost overlap with the analytical solu¬ 
tion except near the boundary r - R. 



FIG. 3. (Color online) The energy shift of the zero-mode solution 
e 0 - For the Neumann BC, it is well explained by the direct formula 
e 0 = -p [Eq. (2.7a)]. For the Dirichlet BC, we find a good fitting 
if we introduce the “effective system radius” R - /3, where the fitting 
parameter ft is determined to be ft = 0.946. 


and i /*2 are separated 15,36 . We consider the system confined in 
a cuboid [~L X , L x \x[-L y , L y ]x[-L z , L z ], and we set aDW per¬ 
pendicular to the je-axis. Henceforth we simply write L x = L. 
The BC can be either Dirichlet or Neumann. Mostly we con¬ 
sider the problem with L y — L, = oo. As shown in Sec. IV, 
strictly speaking, the system with L y — L z — oo has unstable 
modes, i.e., the Bogoliubov equation has the complex eigen¬ 
value. This instability merely reflects the fact that the true 
ground states are the states such that the DW is set parallel 
to the A'-axis, because the surface energy becomes smaller for 
such a configuration. The wavenumbers of unstable modes 
are, however, exponentially small k c ~ s~ aL , and hence we 
can easily suppress these unstable modes by modifying L y , L z 
to be very large but finite sizes satisfying L <sc L yj < -p, which 
makes the wavenumbers of eigenstates discretized and erases 
the unstable modes. 

Let x — d be the position of the DW. By definition \d\ < L 
holds. Let us assume that ipui) occupies the left (right) side 
of the DW, and let pip) be their densities in the uniform re- 


























4 


gion far from both the boundary and the DW. That means, if 
we ignore the detailed profiles near the boundary and the DW, 
the order parameters can be written as \fj\ ~ y[p\9(d - x ) and 
ip 2 ~ \[p20(x - d). When L is large, varying d with fixed p,’s 
corresponds to the smooth sliding of the position of DW with¬ 
out changing the profiles of ip \, ip 2 far from the DW. There¬ 
fore, the differentiation of if /,'s with respect to d with fixed 
Pi’s can be approximated as 


i-d x (x - d) 

(0 (|*-d|»0. 


(2.9) 


with the typical healing length /■'. In particular, if we take the 
limit L —> 00 , we obtain d,i —> —d x . 

The GP equation is given by id, 1/7,■ - - 

( ~ di ~ + 2 J]j=t >2 g, J \Pj\ 1 )Pi, i = 1,2. If L is large, the 

values of p,’s are close to those in the infinite-size system: 
Pi - 2gupi. The Bogoliubov equation can be obtained by 
substituting i/q = i/r, + n,e~ lf7 + v*e 1£ ' to the GP equation and 
linearizing it for («,, v,). 

Now we give our main result on the dispersion relations of 
ripplons in finite-size systems. For simplicity, here we only 
present the result for the case d - 0. The general expres¬ 
sions for d + 0 are available in Sec. IVE [Eqs. (4.77), (4.88) 
with (4.83)]. Let us write the quasiparticle wavefunction as 
(«i, U 2 , Vi, V 2 ) = (i/ifx), u 2 (x), vi(x), v 2 (x))e^ kyy+kzZ) and define 
k = (k 2 +fcr) 112 . Then, the dispersion relation e* and the wave- 
function of the ripplon are given by 


ft = 


2T 0 


tanh kL 


mipi + m 2 p2 


k 2 (k 2 - k 2 ), 


( 2 . 10 ) 


'll \ 


' 772 1 coshk(x + L)lp\ 


(dd4>\\ 

U2 

ft 

—7772 cosh k(x — L) ip 2 

+ tanh kL 

ddih 

Vl 

A2, 

k cosh kL 

-mi cosh k(x + L)ip\ 

, 7722 cosh k(x - L)iA* 

ddP\ 

dd4>%, 

(2.1 


where k c ~ <9(e " L ) is the maximum wavenumber of unstable 
modes mentioned above, and To - f dx( rep¬ 

resents the tension of the DW, recalling the relation Eq. (2.9). 
If we ignore the narrow complex region k < k c , the dispersion 
relation includes the following two cases: 


ft 


27o 


1 mipi + m 2 p2 


VZk 2 

k 312 


(kL « 1) (2.12a) 
(L —> 00 ). (2.12b) 


The behavior e ~ yfLk 2 is consistent with Refs. 18 and 32, 
and the latter case (2.12b) describes the fractional dispersion 
relation 17,18 . The quasiparticle eigenfunction Eq. (2.11) in the 
limit L —> 00 is given by 


Ml' 



1 

( niii//ie kx ' 

111 


d x ip2 

2 T 0 k 

m 2 ip 2 &~ kx 

Vl 


d x ip* 

y 772IP! + m 2 p2 

-mnl/\e kx 

U 2 J 


{d x r 2 ) 


{ m 2 iple~ kx ) 


with recalling dd —> -d x [Eq. (2.9)]. It describes the quasi¬ 
particle wavefunction of ripplons in the infinite system. The 



FIG. 6. (Color online) Dispersion relations of ripplons in the system 
under the Neumann BC. We used the following parameters: 2m 1 = 
2 m 2 = gn = g 22 = y = y = 1, gi 2 = 2.125, and L = 10. 
T 0 is numerically calculated as 2 T 0 - 1.137 with assuming dj = 
-6{6-\x\)d x due to Eq. (2.9). The upper left inset shows the complex¬ 
valued narrow region, and the maximum wavenumber of this region 
is numerically determined as k c — 2.0 x 10~ 6 . The lower-right inset 
shows a plot for larger k’ s. 



FIG. 7. (Color online) Quasiparticle wavefunctions of ripplons in 
the system under the Neumann BC. The parameters are the same as 
those of Fig. 6. The wavenumber is k = 0.2. The ^/-derivative is 
approximated by dj = -6(6 - 1x1)3, due to the relation (2.9). 


former term is the zero-mode solution originated from trans¬ 
lational symmetry breaking. The latter term represents the 
oscillation of relative phases between i/q and i/r 2 and includes 
\[k, indicating that the naive perturbation is impossible. 

Let us see the numerical evidence for the above analytical 
results. We first show the result for the Neumann BC. Fig¬ 
ure 6 shows the numerical verification of dispersion relations. 
An example of quasiparticle wavefunctions is given in Fig. 7. 
The /..-dependence of the quadratic and complex dispersion re¬ 
gions is well illustrated by plotting the /.'-dependence of e*./ k 2 . 
See Fig. 8. 

Our analytical formulas also explain the numerical re¬ 
sults for the Dirichlet BC. As with the case of Kelvin modes, 
we find that the replacement of the effective system length 
L —> L -f3 with p ^ 1.43, and this replacement is used in plot¬ 
ting the analytical formulas. Figure 9 shows the comparison 
of dispersion relations between numerical data and analytical 
formulas with L being replaced by L - ft. Even when we use 
the bare L, a modestly good agreement with the numerical 
data is obtained. However, if we use the modified L - /3, the 
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FIG. 8. (Color online) The log-log plot of k vs e^/k 2 under the Neu¬ 
mann BC. The parameters are the same as those of Fig 6. The three 
solid (dashed) lines show the real (imaginary) part of theoretical for¬ 
mula (2.10) with L = 8, 10,12. The plateau region corresponds to 
quadratic dispersion. The line yJ2To/k corresponds to the fractional 
ripplon dispersion [Eq. (2.12b)]. k c ’s for L = 8, 10 are numerically 
given by k c =* 4.1 x 1(T 5 , 2.Ox 1CT 6 , respectively. k c for L = 12 is too 
small to detect [see Fig. 12 and the paragraph including Eq. (2.14)]. 



FIG. 9. (Color online) The dispersion relation of ripplons in the sys¬ 
tem under the Dirichlet BC with length L = 12. Here, analytical 
formulas [Eqs. (2.10), (2.12a), and (2.12b)] are plotted after replac¬ 
ing L by L — p - 10.57. T 0 is numerically calculated as 2 T 0 - 1.137 
with assuming dj = -0(6 - \x\)d x . The upper left inset shows the 
complex-valued narrow region. k c — 7.2 x 1(T 5 is a numerical fitting 
parameter. The lower-right inset shows a plot for larger fc’s, simply 
showing that the dispersion relation asymptotically comes close to 
that of free particles e = k 2 . 


fitting becomes rather perfect. Figure 10 shows the wavefunc- 
tions of quasiparticle eigenstates. Figure 11 shows the log-log 
plot of q/T 2 , in which the /.-dependence of the quadratic dis¬ 
persion relation becomes visible. The value of /3 is evaluated 
from the plateau region of the data of L - 12 and 16 in this 
figure. 

Here, we give a few additional remarks on the width of 
the complex-valued regions in the dispersion relation, i.e., k c 
in Eq. (2.10). As derived in Appendix B, if we consider the 
system such that 2m\ - 2 m .2 = gn — gn - 1 , g 12 = °°, and 
the average density is given by po = 1 . the /.-dependencies of 
k, for the Dirichlet and the Neumann BCs are given by 


FIG. 10. (Color online) Quasiparticle wavefunctions of ripplons in 
the system under the Dirichlet BC with L = 12. The wavenumber is 
k = 0.2. In using the theoretical formula [Eq. (2.11)], we replace L 
by L - p. The (/-derivative is replaced by dj = -0(6 - \x\)d x . 



FIG. 11. (Color online) The log-log plot of k vs e*/k 2 for dispersion 
relations of ripplons under the Dirichlet BC. The plateau region cor¬ 
responds to the quadratic dispersion. The line y/TToJk corresponds 
to the fractional ripplon dispersion. The results for L = 8, 12, and 
16 are shown. The maximum wavenumbers for the complex region 
are numerically given by k c = 3.9 x 1CT 3 , 7.2 x 10~ 5 , 1.3 x 1CT 6 for 
L = 8, 12,16, respectively. Three solid (dashed) lines represent the 
real (imaginary) part of the analytical formulas [Eq. (2.10)] with L 
being replaced by L -P, P — 1.43. 


Thus, k c in the systems under the Neumann BC decreases 
more rapidly than that under the Dirichlet BC. This relation 
can be also confirmed for finite g 12 with a slight modifica¬ 
tion of the coefficients in exponential factors. See Fig. 12. 
From this figure, we can understand why we cannot find k c 
in the system with the Neumann BC with length L - 12 
in Fig. 8 . We expect k c — 9 x 10 8 from Fig. 12, imply¬ 
ing that the typical eigenenergy of complex-valued region is 
|e[ ~ 0(k 2 .) ~ <9(10 15 ). This is too small to determine k c pre¬ 
cisely in the double-precision calculation. These results are 
consistent with Ref. 18, where the numerical simulations with 
very large U s were performed under the Neumann BC, and 
complex eigenvalues were not found. 


C. Sketch of derivation: summarizing common procedures 


k c oc 


VZe L 

\[Le 2L 


(Dirichlet), 

(Neumann). 


Having presented our main results, we briefly summarize 
(2.14) common procedures of detailed derivations, which will be 
given in Secs. Ill and IV. Even though the mathematical jus- 
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FIG. 12. (Color online) L-dependence of k c , the maximum 
wavenumber of the complex-valued region in the dispersion rela¬ 
tion of ripplons. The physical parameters used are the same as 
other figures. For the Dirichlet BC, the fitting line is given by 
k c = VZexp[-1.05L - 1.89], and that of the Neumann BC is given 
by k c = VTexp[-1.58L - 1.52]. 


this region is independent of the choice of BCs. The behav¬ 
iors of quasiparticle wavefunctions ( u(r ), v(r)) very near the 
boundary L — £ < |r| < L gives no influence to the leading 
order of the dispersion relation e * and is not of our interest 
in the current problem. By the modification (C), the solu¬ 
tion becomes applicable even near a topological defect, i.e., 
0 < |r| < L-(, and thus the effects of zero modes are correctly 
included. For the example of the Kelvin modes, the procedure 

(B) gives the solution F^r) [Eq. (2.4)], and the procedure (C) 
gives Eq. (2.3). For the ripplons, (B) gives cosh k(x ± L) and 

(C) gives Eq. (2.11). See Secs. Ill and IV for detailed deriva¬ 
tions. The evidence of applicability for both Neumann and 
Dirichlet BCs is actually presented in the former part of this 
section. 

Note that, if we consider NGMs concerning spin degree of 
freedom, the terms “density fluctuation” and “phase fluctua¬ 
tion” in the procedure (B) should be replaced by “fluctuation 
of the magnitude of magnetization” and “fluctuation of the an¬ 
gle of magnetization,” respectively. 


tifications for each example are slightly different, the practical 

. ■ . , , ,, III. DETAILED DERIVATION — KELVIN MODES 

procedures are similar. They are summarized as follows: 


(A) First, derive zero-mode solutions having the origin of 
spontaneous symmetry breaking (SSB) in the infinite 
system 32 . 

(B) In the intermediate region far from both topological de¬ 
fects and the boundary, where the asymptotic form of the 
order parameter becomes almost exact, derive the finite- 
wavenumber solution of the Bogoliubov equation. In 
such a region where the local structure of the order pa¬ 
rameter is ignorable, the density fluctuation (~ u + v) be¬ 
comes irrelevant compared to the phase fluctuation (~ u - 
v), and hence the differential equation becomes solvable. 
Here, the integration constants are fixed by assuming the 
Neumann BC lim r ^boundary n • Vu(r) = n ■ Vv(r) = 0. 

(C) Make a minimal modification to the solution obtained in 
(B) to include the exact zero-mode solutions derived in 
(A) to take into account the local structure near the topo¬ 
logical defects. 

(D) Using the solution constructed in the above way, calcu¬ 
late an eigenenergy tk solving the Bogoliubov equation 
by using the techniques in Ref. 32. 

Here, we emphasize that the use of the Neumann BC in the 
procedure (B) does not mean that our result is not applica¬ 
ble for other BCs, e.g., the Dirichlet BC. The purpose of (B) 
is to obtain the quasiparticle wavefunctions in the asymptotic 
region where the behavior of the order parameter becomes al¬ 
most uniform. Since u, v are linearized fields of the order pa¬ 
rameter, they also should obey the same uniform boundary 
condition, and hence the Neumann BC is most suitable for 
this purpose. To be more concrete, let |r| be a distance from a 
topological defect and let £ and L he a typical healing length 
and the distance between the defect and the boundary, respec¬ 
tively. Then, the solution obtained in (B) is quite applicable in 
the intermediate region £ <§: |r| < L - £, and the behavior in 


Thus far, we have presented our main analytical formulas 
and their numerical verifications for Kelvin and ripple modes. 
In this and next section, we provide the complete derivations 
of these formulas. 


A. Fundamental equations and zero modes 


The energy functional of the one-component BEC in the 
dimensionless form is given by 


H-fiN = J d 3 r- (|V<A| 2 + |<A| 4 - 2|^| 2 ) . 

The stationary GP and Bogoliubov equations are 

(-V 2 - 2 + 2|^| 2 )«A = 0, 
/-V 2 - 2 + 4\i//\ 2 2i// 2 \ tu\ _ lu\ 

\ V 2 + 2 - 4|i/r| 2 j yvj ~~ U7 ' 


(3.1) 


(3-2) 

(3.3) 


Since we are interested in the vortex solution with the vortex 
charge n = 1, we set if/ - where f(r ) is a non-negative 

function having the asymptotic form /(oo) = 1. Then the GP 
equation becomes 

-/"-- + 4-2/(l-/ 2 ) = 0. (3.4) 

r r- 


Henceforth we write the solution in the infinite-size system as 
f<x>(r). The asymptotic solution is given by 


= 1 - 4? 


9 

32 ^ 


+ 0(r 6 ). 


(3.5) 


The expansion at r — 0 can be also obtained, given by 


foo(r) = ar- -r 3 + 


a + 4t/ 3 


r 5 + 0(r 7 ), (3.6) 
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where a ^ 0.82 is a constant determined numerically. 

In the infinite-size system, the GP equation has a symmetry 
such that “il/(x,y,z) is a solution” «-> + xo,y + yo,j ) el9 

is also a solution”. Differentiating the GP equation by 9, xq, 
and yo, we obtain the following SSB-originated zero mode 
solutions 32 for the Bogoliubov equation: 


which holds for any “Bogoliubov-hermitian” operator 32 , and 
can be regarded as an analog of self-adjointness for hermi- 
tian operators. Using these inner products and analog of self¬ 
adjointness, we can construct a perturbation theory in a similar 
way to that of ordinary hermitian operators 32 . 


Wphase 



Wx-trans 


I d x if/\ 

\d x i!/*J ’ Wv ' trans 



(3.7) 


B. Type-II dispersion coefficient in finite systems 


As shown in Ref. 32, w p h as e is cr-orthogonal to the other two 
zero modes, so it solely yields a type-I NGM, which is the 
Bogoliubov phonon. On the other hand, wvt ram and vv v .t rans 
are not cr-orthogonal and becoming a pair yielding one type- 
II NGM, the Kelvin mode. We can construct a positive-norm 
zero-mode solution becoming a seed of type-II mode by their 
linear combination, which is given by 32 


Wq - W^.trans _ Wy-trans 


f + h. 

J OO r 

e- 2iS (/; - £). 


(3.8) 


Then wo becomes the seed of the positive dispersion branch. 
The same solution was also derived by Pitaevskii 5 . vv A . t rans + 
inytrans has negative norm and yields the negative dispersion 
branch. 

The Bogoliubov equation can be decoupled for different an¬ 
gular momenta by setting (u, v) = (n(r, z)e lfl , v(r, z)e~ lfl )e lm61 , 
m 6 Z. We are further interested in the solution propagating 
in the "-direction. So we set (n(r, z), v(r, z)) = e*** z (M(r), v(r)). 
The resultant equation is 


Henceforth we consider the finite-size systems. Let the sys¬ 
tem be an infinitely-long cylinder with finite radius R. The 
BC at r — R is arbitrary and does not give an influence to the 
following argument. In a finite-size system, the translational 
symmetry no longer exists and hence w. v _ t rans and wy tr ans do not 
become the exact zero-mode solutions. Let us see how these 
zero-mode solutions are modified in finite-size systems. 

We solve the Bogoliubov equation using the expansion w.r.t 
the parameter a R 1 . Then a - 0 corresponds to the 
infinite-size system R = oo and finite a corresponds to finite- 
size systems. Let us write £ = r/R = ra. Then £ can 
take a value in the closed interval [0,1], Let us further write 
fig, a) f(£/a). We henceforth use the prime symbol to 
express the ^-derivative, e.g., f' - %■ Then the GP equation 
(3.4) becomes 

W -/"-j + i)-2f0-p) = 0. (3.17) 


6 (") = (//„ + rf)(j) (3.9) 

with k = \k z \, cr = diag(l, -1) and 

t 1 (in + l) 2 t 

[ffoht =-d;--d r + -—— 2 + 4/ 2 , (3.10) 

r r L 

[H 0 ] l2 = -[H 0 \ 2l =2f 2 , (3.11) 

, 1 (m — l) 2 t 

[Hohi = d; + -d r - y— + 2 - 4/ 2 . (3.12) 

r r- 

The Kelvin mode with positive dispersion exists in the sec¬ 
tor m = -1, since it contains the zero-mode wq [Eq. (3.8)]. 
Henceforth we consider only this sector. The asymptotic be¬ 
havior of zero-mode solution is given by 


Let us seek a solution in the form of n-expansion: f - fo + 
& 2 fi + < 2 4 / 4 + • • • ■ Note that the expansion around a - 0 is 
rather sensitive and only meaningful in 0 < £ < 1. At £ = 0 
and 1 , the expansion is pathological and we do not consider 
it. Here, we are only interested in the intermediate regions far 
from both the vortex and the boundary. The GP equations for 
each order then become 

ff°: /o(l~/o) = 0, (3.18) 

: - fo - T + t - 2/ 2 (i - 3/;-) - 0, (3.19) 

£ r 

- 4 + f - 2/4(1 - 3/^) + 6 /o/? = 0. (3.20) 


U = f£,(r) + 


foc(r) 



v = fL(r) - 


foo(r) 



(3.13) 

(3.14) 


The cr-inner products between two quasiparticle wavefunc- 
tions Wj = ( Ui(r ), v,(r)), i = 1,2 is defined by 


Oi,w 2 )<r := f rdr(u\u 2 - vjv 2 ). (3.15) 
Jo 


Here we omit the 0-integration, which merely gives the factor 
2n in this problem. Hq satisfy the following property 


(w\,HqW2)o- - (HqWi, W2)cr, (3.16) 


The solution satisfying the asymptotic condition f(r —> oo) = 
1 is given by fo — 1 , and fa, fa,... are determined iteratively: 

/o=l, /! = -4j. A = - 3 |j. - < 3 2 » 

Thus we have 

- 2 -^ 4 ? +0{a6) - (3 - 22) 

This is just the revisit of Eq. (3.5). 

Next we solve the Bogoliubov equation by the same expan¬ 
sion. The Bogoliubov equation rewritten by £ and a is given 








by 


a 


2 



a \v + 


£ 


j — 2(1 — 2 / 2 )m + k^u + 2/ 2 v = eu, (3.23) 
■ + 2(1 - 2 f)v + k 2 v - 2fu = ev. (3.24) 
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Substituting Eqs. (3.30) and (3.32) to Eq. (3.26) and compar¬ 
ing the coefficient of £ l in both sides, we obtain eo ,2 = -2. 
Thus, the energy shift of the zero-mode solution in the finite- 
size system becomes 

eo = + 0(R 4 ) (Neumann BC), (3.33) 


Here we again note that the prime represents the differentia¬ 
tion by £. 

We first consider the zero-wavenumber case k - 0 and ex¬ 
amine the energy shift of the zero-mode solution wq due to the 
finite-size effect. Let eo be the energy shift of the zero-mode 
solution, and let us expand it as eo = eo.o + a 2 60,2 + o 4 e 0,4 + • ■ ■. 
We already know that the eigenvalue of vvq in the infinite sys¬ 
tem (a = 0) is zero: eo.o = 0. We also expand the quasiparticle 
wavefunctions in the same way: 11 = uq + a 2 U 2 + ■■■ , v = 
vo + a 2 v ’2 + • • ■. The zeroth- and second-order equations are 
then given by 


M 0 + V 0 = 0, 

(3.25) 

) M 0 

-—r + 2 (M 2 + V 2 ) = 60 . 2 M 0 , 

£ 

(3.26) 

3v’o 

- ~pr ~ 2(m 2 + v 2 ) = e 0 , 2 V 0 . 

£ 

(3.27) 


as with Ref. 31. If we use the Dirichlet BC, we find a little 
larger correction due to the boundary effect: 

e 0 = - J- + 0(7T 3 ) (Dirichlet BC), (3.34) 

though the leading order is the same. 

The solution (3.31) well describes the numerical solution in 
the region 0 « r < R, but it diverges at r — 0. This artificial 
divergence is caused by the fact that the a-expansion is valid 
only for £ 6 (0,1). In order to get the correct behavior near 
the vortex core r = 0, we heuristically replace the divergent 
term r~ l by the zero-mode solution of infinite systems, i.e., 
Eqs. (3.13) and (3.14). This replacement is good if the system 
size R is sufficiently large, because the profile of quasiparticle 
wavefunctions near the vortex core is almost the same with 
those of infinite-size systems. Thus, we obtain the modified 
zero-mode solution in the finite-size system as 


Thus we obtain mo+vo = 0, which justifies ignoring the density 
fluctuation in the procedure (B) of Sec. IIC. Taking the sum 
and difference of Eqs. (3.26) and (3.27), and using vo = -uq, 
we obtain 


~ J + p = °- (3-28) 

2 un 

—-pr + 2(W2 + V 2 ) = eo, 2 «o- (3.29) 

£~ 

The solution of Eq. (3.28) is given by u = c\£+cn £~ 1 . Follow¬ 
ing the procedure (B), we fix the coefficient by the Neumann 
BC u' 0 (£ -> 1) = 0. Thus, 


uq--vq-£+~. (3.30) 

If we go back to the original variables, r and R, this solution 
can be rewritten as 


«o - -v o - - + Tpr- (3.31) 

r R~ 

While the term 2 corresponds to the expansion of the zero¬ 
mode solution in the infinite-system Eqs. (3.13) and (3.14), 
the latter term jp exists purely by the finite-size effect. This 
term is necessary to obtain the energy shift eo, 2 - 

Let us find the expansion coefficient eo, 2 - To derive this, 
we focus on Eq. (3.29) in the region 0 < £ «: 1. Using the 
next leading orders of Eqs. (3.13) and (3.14), in the region 
0 <f«l, the leading order terms of M 2 and V 2 are given by 

112 = ^3 + V 2 = ^3 + 0 ^"‘)- 


Wq 



JL + f + w) 
V 2 I 


(3.35) 


where the factor 4= is a normalization factor. This modifica- 

V2 

tion corresponds to the procedure (C) in Sec. IIC. 

Using Eq. (3.35), we can calculate the coefficient of type-II 
dispersion. Let us solve the Bogoliubov equation perturba- 
tively: 


(Hq + crk 2 )(wo + k 2 W 2 H-) 

= (e 0 + k 2 ei + ■ ■■ )(w 0 + k 2 w 2 + ■■•)■ (3.36) 


The zeroth and the second order equations are 

w 0 = e 0 w 0 , (3.37) 

HqW 2 + CTWo = eoW 2 + 62 Wq. (3.38) 


Here we already know eo = — ^ + 0(R 4 ). Note that 62 in this 
^-expansion is different from eo ,2 appearing in the cr-expansion 
of eo- Taking the cr-inner product between wo and the second 
order equation, and using (3.16), we have 


(wq,otvo) 0 . = Ip R rdr ( ll l + vg) 
( w 0 , w 0 )fr J 0 R rdr(u 2 -v 2 ) 


The denominator is evaluated as 

-R 


(wo, Wo)c 


-f 


dr 


2fL(r)\U(r)+p 


log R 


R 2 


(3.39) 


(3.32) 


= 1 + 0 


(3.40) 
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where the orders of each term are evaluated using Eq. (3.5): 

f dr[2/j>)/ ra (r)] = U(R) 

Jo 

X R r*R 

dr\fL(r)r 2 \ ~ J dr 


2 = 1 + <9(7C 2 ), (3.41) 


log R. (3.42) 


In both cases, however, the leading term is the same. 

As shown in Eqs. (3.34) and (3.52), the Dirichlet BC gives 
a little larger deviation from the leading order term compared 
to the Neumann BC. As discussed in Sec. II, these deviations 
are well included by the effective replacement 

R^R-P, p = 0.946, (3.53) 


Thus, wo is normalized up to 0(R 2 log R) terms. The numer¬ 
ator is given by 


(wo, CTW 0 )o- 



2r/oo(r) 

R 2 


+ rfL(r) 2 + 


foo(r ) 2 


(3.43) 


The leading orders of each term are given by 


f 


dr 


r 

¥ 


2 rfoo(r) 


1 

— ¥)■ 


f 

X R poo 

d r[rfL(r) 2 } = d r[rf^(r) 2 ] 


(3.44) 

(3.45) 

+ 0(R\ (3.46) 


and 


f 


dr 


foc(r)~ 


= [U(r ) 2 log -f Q dr \ 2 Ur)f'Jr) log r\ 

dr[2./„(r)/l(r)]ogr] 
o 


O 


log R\ 

r 2 r 


(3.47) 


Here, /oo(r ) 2 logr| r= o vanishes since foo(r) — ar [Eq. (3.6)]. 
Thus, we obtain 


(wo, crwo)cr = log R + -+ rj + O 


log R\ 
R 2 )’ 


(3.48) 


X oo 

dr [rfLirj 2 - 2Mr)f„(r) log r] - 0.227. (3.49) 


A closed form for this ;/ is not known. 
Summarizing, we obtain 


R 2 


+ Ak 2 + 0(k 4 ) 


(3.50) 


with 


A = log R + — + rj + O 


log R 
R 2 


(Neumann BC). (3.51) 


If we use the Dirichlet BC [f(R) = u(R) = v(R) = 0], the 
profile of quasiparticle wavefunctions near the boundary r - 
R deviates from wo = (»o, Vq ) 7 - This deviation yields a little 
larger correction: 


A = log R +^+n + 0(R X ) (Dirichlet BC). (3.52) 


where the value of p is determined by numerical fitting of eo 
(Fig. 3). The physical meaning of this replacement is as fol¬ 
lows. Since the order parameter is suppressed near the bound¬ 
ary, the effective radius of the system becomes about a healing 
length shorter than that of the Neumann BC. See also Fig. 5. 

The formula obtained here explains the numerical results 
very well for small wavenumbers in finite-size systems. How¬ 
ever, we cannot take the limit R —» oo in this expression. In 
the next subsection, we derive an interpolating formula valid 
even for R - oo. 


C. Interpolating formula, derivation of e - k 2 log k 

Now we consider the finite-wavenumber case of Eqs. (3.23) 
and (3.24). Since we are interested in the region such that 
kR ~ (9(1), we expand the wavenumber as k - ka + ■ ■ ■. 
The energy and quasiparticle wavefunctions are expanded in 
the same way with the previous subsection: e = eg + a 2 e 2 + 
■ ■ ■ , (m, v) = (no, vo) + a 2 (u 2 , V 2 ) + • ■ •. Then, the zeroth-order 
equations are 


2 (« 0 + v 0 ) = eo«o, (3.54) 

-2 (« 0 + v 0 ) = e 0 v 0 . (3.55) 

In order for these equations to have a nonvanishing solution, 
det( 2 + 2 ° 2 -f (: ) = 0 is necessary. Thus we obtain eo = 0 and 
«o + vo = 0 , which again gives the justification for the pro¬ 
cedure (B) in Sec. IIC. The second order equations are given 
by 


,, Uo 

-M 0 - — - -7 + £ n 0 + 2(n 2 + v 2 ) = e 2 n 0 , (3.56) 

? r 

v o + -7 - ^pr ~ k 2 v 0 - 2(n 2 + v 2 ) = e 2 v 0 . (3.57) 

k k 

Taking the sum of these two equations and using no + vo = 0, 
we obtain 


~ + ? + pM ° =0> (3 - 58) 

which is just the modified Bessel differential equation. Thus 
the solution is given by no = c\I\(kp) + (kp). Again, 

following the procedure (B), imposing the Neumann BC 
lim^i u'{p) = 0 , we obtain 

u 0 = ~k[K 1 (kt)+xW 1 (kt)\, (3.59) 

, r . Ko(k) + K 2 (k) 

Io(k) + I 2 (k) 


(3.60) 
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If we go back to the original variables r and R. this solution 
can be rewritten as 


«o = F k (r) := k[K\(kr) + X (kR)h(kr )\. (3.61) 


This F k (r) has a few notable properties. If R + oo, it has a 
Taylor series around k = 0: 


N 1 r k 2 r 

FjAr) — — + —— +- 

ky ’ r R 1 8 


-7 + + 4 log — I + 0(k), (3.62) 


which implies that the naive perturbation works out well if the 
system size is finite. On the other hand, if R = oo, the function 
X (k) has the asymptotic behavior 


r k) _|p - T-f- log § + 0 (k 2 ) (~k« 1) 

\ne- 2k [l + ± + 0(k- 2 )] (k» 1). 


(3.63) 


Hence, Hindoo F k (r) = kK\ (kr). which does not have a Tay¬ 
lor series since K\(kr) includes the logarithmic term. Thus, 
we cannot use the naive perturbation theory in the infinite- 
size system. We mention that the solution u - K\(kr) was 
also found by Pitaevskii 5 . 

Now, following the same procedure with the previous sub¬ 
section, we modify this solution in order to avoid the artificial 
divergence at r = 0. Namely, we use the following modified 
quasiparticle wavefunction: 


w k := 


/ Uk(r)\ 
\Vk(r)J 


V2 


Fk(r) - ± H 

i-T,(r) + i-(^ 


fLirj) ' 


(3.64) 


This expression just gives Eq. (2.3) up to a factor. If we set 
k - 0 in this expression, we again obtain Eq. (3.35). 

Let us calculate the eigenenergy e k of w k by solving the 
Bogoliubov equation [Eq. (3.9)] 


(H 0 + crk 2 )w k = e k w k . (3.65) 

Taking the cr- inner product between this equation and wq, we 
obtain 


where 


h = 
h 


|| d+ 

= f rdr,C(r)[F k (r) - F 0 (r)], 
Jo 

r R 

h = I rdr [u 0 (r) 2 + v 0 (r) 2 ], 

*-/* 

k-f 


rdrF 0 (r)[F k (r) - F 0 (r )], 


dr[f«,(r)-l][F k (r)-F 0 (r)]. 


(3.69) 

(3.70) 

(3.71) 

(3.72) 

(3.73) 


The integrals 7i and It, are ^-independent and already evalu¬ 
ated in the previous subsection [Eqs. (3.40) and (3.48)]: 


7, = 1 + O 


log R\ 

r 2 r 


It, = log R + ^ + Tj + O 


log R 
R 2 


(3.74) 

(3.75) 


If we perform the order evaluation by regarding k = 0(R 1 ), 
h and R are shown to be ignorable: 

^(logft) 2 \ / (log R) 2 ^ 


h-0 


R 2 


h - O 


R 2 


(3.76) 


7 4 can be symbolically integrated as 


h = 


X(kR) I I 0 (kr) + 


(log r + K 0 (kr )) • 


2 h(kr)\ 

R 2 ) 

r-K 2 (kr) r 2 r 

R 2 


.4 1* 

R 2 4^J 0 


2 , kR 5 

-- A -y-xW), 


(3.77) 


where the behaviors Ko(kr) + log r = —y - + 0(r 2 ) and 

K 2 (kr) - + 0 ( 1 ) are used. 

Summarizing, the dispersion relation of the Kelvin mode is 
given by 


(3.78) 


e k = k 2 1 - log - + 77 - y - x(kR) 


Q — 60 + 


k 2 


- eo + k 2 


(wo, crwk)^ 

(wo, W k )o- 

r R 

J 0 rdr [u 0 (r)u k (r) + v 0 (r)v jt (r)] 

pR 

J 0 rdr [u 0 (r)u k (r) - v 0 (r)v 4 .(r)] 


(3.66) 


This formula includes the following two limiting cases: 


e k - 


_ ^2 + ^ ^log R + — + TjJ (kR « 1) (3.79a) 
k 2 1- log ^ - yj (R -> 00 ). (3.79b) 


We already know eo - -^r [Eq. (3.33)]. Let us calculate the 
inner products. We write 

rdr [u 0 (r)u k (r ) - v 0 (r)v*(r)] = 7] + I 2 , (3.67) 

rdr [uo(r)u k (r) + \’o(r)v k (r)] = 7 3 + 7 4 + I 5 , (3.68) 



The case kR <s 1 revisits Eq. (3.50) with (3.51). The latter 

case gives the non-integer dispersion e - k 2 log k. which was 

first shown in Ref. 5. Taking the limit R —» 00 in Eq. (3.64), 
the quasiparticle wavefunction of Kelvin modes in the infinite 
system becomes 


[u k (r)\ = _L ( kK ^ - 7 + + f~(r) ) 

V’k(r)) V2 \-kKJkr) + 7 - + fcoir)} ' 


(3.80) 
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Finally, we would like to give a few perspectives on the 
higher-order corrections of the dispersion relation Eq. (3.78) 
and its infinite limit Eq. (3.79b). In deriving this formula, we 
have ignored the terms I 2 and Is by assuming k = 0 (R '). 
Since they vanish at k = 0, this ignoring is not bad even in 
the infinite system R = 00 , if k is small. Indeed, the nu¬ 
merical result with R = 100 given in Fig. 1 shows that the 
formula (3.79b) is good for 0 < k < 0.1. However, if we 
are interested in the next leading order term of the formula 
(3.78), we must include contributions from I 2 and Is- The 
emergence of 0((\og R ) 2 / R 2 ) terms implies that, if these terms 
are treated with mathematical care, they will become of or¬ 
der 0 (k 2 (\ogk) 2 ), meaning that the next leading term of the 
dispersion relation Eq. (3.79b) would be given by A: 4 (log k) 1 . 
However, at this time, we do not have a derivation for this con¬ 
jecture and a possible finite-size generalization. This is left to 
be an open problem. 


IV. DETAILED DERIVATION — RIPPLONS 


In this section we provide the detailed derivations of ana¬ 
lytical formulas for ripplons presented in Sec. II. 


A. Fundamental equations and ground states in ID systems 

We first consider the ground state of the one-dimensional 
system with length 2 L 


J-l \ 2 m 1 2 m 2 


+ glll(^ll +2gl2|<Al| Itfol + g22\^l\‘ 


(4.1) 


with fixed particle numbers N, = f dxl^j 2 , i = 1,2. Though 
the result of this problem is well-known 15,36 , we review it 
in order to introduce the variable d [Eq. (4.10)], having the 
meaning of the position of the DW. The discussion given be¬ 
low holds regardless of whether the BC at x — ±L is of Dirich- 
let or Neumann. 

Let us assume that the system length 2 L is sufficiently large 
compared to the typical healing length of the order param¬ 
eters and hence the energies of bulk condensates are much 
larger than those of surfaces and boundaries. (We can intro¬ 
duce four kinds of healing lengths in this system as seen in 
Appendix A.) Assume that two condensates \j/\ and ifj 2 are 
separated, and i/qp) occupies the left (right) side of the box 
with length ii( 2 ), where L\ + Li - 2L. Then, the energy of 
this state is given by 


N 2 N 2 

H separated — gll~ ^ g22 2T: • (4-2) 

L,\ AL, — L,\ 


Minimization of ^separated with respect to L\ yields 


U_ = . = 1 2 

2 .L sJgnNi + '\fg22N2 
tt _ + yfg 22 ^ 2) 2 

H separated • 


(4.3) 

(4.4) 


On the other hand, as another ansatz, the energy of the uniform 
mixture of t/q and i/q is given by 


g\\N 2 + g 22 ^\ + 2 gnN\N 2 
Hmixed - -^-, (4.5) 

which does not have an additional parameter to be optimized. 
The energy difference between these two states is given by 

TT TT _ NlN 2 (gl2 - '\JgWg22) fA ^ 

^mixed ^separated — ^ 

Thus, if gi 2 > y[gng 22 , the ground state is given by the state 
such that i/q and <p 2 are separated. 

Henceforth we only consider the separated case. The den¬ 
sities of these condensates are given by 



Vc?ll^l + ^/g 22^2 
2 Lx/gJi 


i = 1 , 2 . 


If we introduce 


(4.7) 


P ■= 



the densities can be rewritten as 


(4.8) 


P = \[guP\ = \fg 22 P 2 - (4.9) 


This relation also holds in the infinite-size system due to the 
momentum conservation law (see Appendix A). The position 
of the DW is given by 


L( yfgnN\ - yfgziNi) 
y/gnNi + VS 22 AV 2 


(4.10) 


We can use p and d as system parameters instead of N\ and 
N 2 - The relation between them are 


Ah 


P(L + d) ^ _ p(L - d) 
yfgu ’ “ y/g 22 


(4.11) 


Henceforth we regard i/r,’s as functions of these parameters 
instead of N t and N 2 , that is, they are considered as a function 
tpi - il/j(x,p,d). If the system length 2 L is sufficiently large 
and the DW is located far from the boundary (i.e., |d ± L\ is 
much larger than the typical healing length), changing d with 
fixed p implies a smooth sliding of the DW almost without 
changing the profiles of i/q, t/q far from the DW. If gn = g 22 , 
the story becomes a little simpler; since p oc N\ + N 2 and 
d cc N\ - Ah, the sliding of the DW occurs by changing the 
imbalance of the particle numbers N\ - AA 2 with fixing the total 
number N\ + AA 2 . In the general case gn + g 22 , however, fixing 
p does not mean fixing the total particle number. 
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From the above physical interpretation, the differentiation 
with respect to d with fixed p is approximately given by 


{x " d) 

dd |0 (\x -d\» £), 


(4.12) 


with f being the typical healing length. In particular, if we 
take the infinite-size limit, we have 


lim tt; 

£->°° dd 


d_ 

dx 


(4.13) 


B. SSB-originated zero mode solutions and overview of 
calculation 


Now let us consider a three-dimensional system. We con¬ 
sider the system such that the length with respect to the x- 
direction is 2 L and those with respect to the y- and "-directions 
are infinite: L y - L- - oo. As shown below, if L y — L z — oo, 
the Bogoliubov equation has the complex eigenvalue. In 
other words, the system has unstable modes. However, the 
wavenumbers of unstable modes are shown to be exponen¬ 
tially small k ~ e aL , therefore we can easily eliminate 
these unstable modes through discretization of wavenumbers, 
which is realized by modifying L y , L z to be very large but fi¬ 
nite sizes. 

Let us consider the GP and Bogoliubov equations. Assum¬ 
ing translationally-invariant configurations along the y- and in¬ 
directions, the GP equation is reduced to 

{-Vi-2tr l d 2 x + 2g u \ti\ 2 + 2gn\il'2\ 2 )ti=0, (4.14) 

(~H2 - ^d 2 x + 2g2ll<Al| 2 + 2g22\tl'2\ 2 )l/'2 = o. (4.15) 


For the DW solution, if/\ and 4*2 can be taken as real-valued 
functions up to overall phase factors. The chemical potentials 
are the functions of system parameters: p- t — pdp, d), and they 
are determined via the condition f dx|iA 1 j 2 = N,-. If L is large, 
they are almost the same with those of the infinite-size sys¬ 
tem: pi - 2gupi = 2 /? (see Appendix A). Therefore, the 

^-dependence of //,’s is expected to be very small for large L. 

For the Bogoliubov equation, assuming the plane-wave 
solution in the y- and z- directions, we set Uj(x,y,z ) = 
Ui(x)e' l - kyy+kzZ \ Vi(x,y,z ) = Vi(x)e' < ' k > y+kr - z \ i = 1,2. We then 
obtain 


(H 0 + Mok 2 ) 


'lt\ 


'u\ 

U2 

= € 

M2 

Vi 


Vl 

4’2, 


yi, 


(4.16) 


where k = 


-yjtf + k\. Mo = diag( 


1 1 -1 

2 m 1 ’ 2m2 ’ 2 m 1 ’ 


-1 

2/722 


), and 


H 0 


I F 0 G 0 \ 

l-c; -f;J 


(4.17) 


with 


F 0 


Go 


2mi 2 m 2 ^ 2 ) 


4gn|(Atl 2 + 2g l2 \4/ 2 \ 2 

^g\l4*\4*2 


Zgnil'dl'l 

4g22|lA2| 2 + 2 g 


12 


l<Ai|- 


(4.18) 


/ 2gii(A 2 2g l2 i// l i// 2 \ 
\2gi 2 i//i(//2 2g 22 i//l ) 


(4.19) 


Note that the kinetic energy term is not cr = diag( 1,1,-1, — 1), 
because the masses are generally different: m\ 4- m 2 . The cr- 
inner product between two quasiparticle wavefunctions w,- = 
(m,-i, uq, v ;1 , Va) T , i — 1,2 is defined by 

(Wl ,W 2 )<r = J' dx(u* n lt 2 l + U n ll 22 - Vj[V21 - Vi 2 V22). 

(4.20) 


Hq and Mq satisfy the “Bogoliubov-hermitian” property 32 : 

(x, Hoy) a = (H 0 x, y)o-, (4.21) 

(x, M 0 y)cr = (Mqx, y) lT . (4.22) 

Let us discuss SSB-originated zero-mode solutions 32 . In 
the infinite-size system, if (>fj\(x,y,z),4/ 2 (x,y,z)) is a solution 
of the GP equation, (e lfl| if/\(x + xq, y,z), e 102 i ^ 2 (x + xo,y, z)) is 
also a solution. By differentiating the GP equation with re¬ 
spect to 0i, 0 2 and xo, we have the following zero-mode solu¬ 
tions: 



r 


' 0 ' 


yr 


0 


4*2 

_ d 

4*2 

W\ = 

-r x 

, W2 = 

0 

•> Wtrans — ~7T~ 

OX 

4 *; 


{0 J 


W 2 ) 


[ r 2 ) 


However, if we consider a finite-size system, only w\ and vv 2 
are exact zero-mode solutions and w tr ans is no longer a solu¬ 
tion since the translational symmetry is absent. In the finite- 
size system, the generalized eigenvector Zd, derived in the next 
subsection, plays an alternative role to Wtrans- 

Since these two modes are cr-orthogonal to each other 
(vvi,W 2 )o- = 0, we conclude that the system has two type-I 
NGMs and no type-II NGM appears by following the gen¬ 
eral theory constructed in Ref. 32. At first glance, this fact 
would seem contradictory to the fact that the ripplon has a 
type-II dispersion in a finite-size system 18 - 32 . This apparent 
paradox can be resolved in the following way: the gapless 
mode corresponding to the ripplon indeed has a linear disper¬ 
sion e - ak in finite-size systems. However, the coefficient a 
is an exponentially small complex number. If we ignore this 
exponentially small region k < 0(e / /2 ), the dispersion rela¬ 
tion for k < 0(L l ) is well described by e ~ x[Lk 2 , as shown 
in Refs. 18 and 32. Furthermore, if k becomes a little larger, 
the dispersion relation becomes e ~ k 2 ^ 2 . These three differ¬ 
ent behaviors in different wavenumber scales will be solely 
explained by one formula in Eqs. (4.77) and (4.83), which 
are the goal of this section. Henceforth, we solve the Bogoli¬ 
ubov equation in the three ways shown in Table I to derive 
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TABLE I. A list of approximations and derivable dispersion relations 
for the ripplons in finite-size systems. The naive perturbation, two- 
state approximation, and /c-dependent two-state approximation are 
discussed in Subsecs. IV C, IV D, and IV E, respectively. 



e oc i ck 

e oc yfZk 2 

e oc k 3/2 

Naive perturbation 

y 



Two-state approximation 

y 

y 


k-dependent two-state approximation 

y 

y 

y 


the above-mentioned three behaviors. Even though the last 
method provided in Subsec. IV E gives the most general and 
important result, the former methods treated in Subsecs. IV C 
and IV D are necessary to formulate the last method. So we 
need all three formulations. 


C. Naive perturbation — type-I complex dispersion 

We first solve the Bogoliubov equation by a naive perturba¬ 
tion theory and find the complex-coefficient type-I dispersion. 

Since w 1 and W 2 are the seeds of type-I NGMs, there must 
exist generalized eigenvectors satisfying Hqh w ; according 
to Ref. 32. Such vectors can be found by differentiating the 
GP equation with respect to the system parameters 32,38 . The 
differentiation with respect to p and d yields 

H 0 Zp = pipWi + p2 P W2, (4.24) 

Hold - P\dW\ + P2dW2, (4.25) 

where 


*1 





dz 

dz 

(4.26) 

Vx 

r 2 ) 

, z P :-—, Zd 

dd’ 

dpt dpi 

dp ’ dd' ' 

= 1,2. 

(4.27) 


Let us define the following notation for later convenience 


This implies that ptd vanishes if we only take the leading order. 
A rigorous evaluation of pu,p 2 d is not easy, but the typical 
behavior is given by 

p id ~ Le'^, p ul < 0, p 2d > 0, (4.32) 


where is the typical healing length of order parameters and 
a is an 0(1) constant. For the special case gn = +°°, we 
can rigorously derive the behavior in Eq. (4.32), because the 
two condensates are completely separated and hence the GP 
equation reduces to that of a single-component BEC. See Ap¬ 
pendix B. We can also find similar behaviors for finite gn 
from numerics. As we see below, these small ptd ’s cause a 
very narrow complex eigenvalue region in the dispersion rela¬ 
tion. Since ptd’s are very small, we often ignore higher-order 
terms of p,d’s in the following calculation. 

Because of Eqs. (4.25) and (4.32), the generalized eigen¬ 
vector Zd is an “almost” zero-mode solution if L is large. 
In particular, using Eq. (4.13), it exactly reduces to the zero 
mode solution due to the translational symmetry breaking in 
the infinite-size limit: 


Zd 


d 

—Wtrans — 2T" 

OX 


yr 

1 A 2 

l r 2 ) 


{L 


00). 


(4.33) 


This relation implies that Zd plays an alternative role to w tr ans 
in finite-size systems. 

Let us derive the eigenvectors and eigenvalues of the Bo¬ 
goliubov equations (4.16) by solving it perturbatively 32 . Let 
us look for the eigenvector and eigenvalue by the expansion 

f = fo + fi* + f2* 2 + -”, (4.34) 

e = eik + e 2 k 2 + ■■■ (4.35) 


with 


£0 = Q\W\ + CI 2 W 2 , (4.36) 

ft = b\z\ + b 2 Z 2 - (4.37) 

The zeroth order equation Hq£q = 0 holds identically. From 
the first-order equation = e\£a, we obtain 

b 1 = eiai, b 2 = e\a.2. (4.38) 


[ ’ pd — p ~dd ~ — p ~dd' 
Then, if we introduce 

_ [Z,pi\pd ^ _ ldi,z] p d 

]pd’ " [J^\,^^2]pd , 

they satisfy 

HoZi = Wt, i - 1,2. 


(4.28) 


(4.29) 


(4.30) 


As already mentioned, if the system length L is sufficiently 
large, pt can be approximated by those of infinite-size sys¬ 
tems: pt =: 2 gupi ^ 2yfgitp. (See Appendix A.) Therefore, 

(4.31) 


The second-order equation is given by + T/ 0^2 = e 2 ^o + 
ei<(i. Taking the cr-inner product between this equation and w ( - 
gives 

where W and G are 2x2 matrices whose components are 
defined by 


[W]tj = (Wj, MqW j)o-, [G]tj = C Wi, Zj) a ■■ 


They can be calculated as 


W = 

G = 


(Ni/mi 0 \ 

( 0 N2/ni2)’ 

1 t[Nl,P2\pd 
[pi ,P2]pd \[7V 2 ,jU2 ]pd 


[Mu ^Vi]pd\ 
[fJl,N2]pdj 


(4.40) 


(4.41) 


Pip - 2 yfgii, ptd - 0, i - 1,2. 


(4.42) 
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Note that the entries of G can be also written as [G],q = 
(HoZi,Zj)o- by Eq. (4.30), and hence G is hermitian due to 
Eq. (4.21). Furthermore, Eq. (4.42) shows that G is real, hence 
G is a real-symmetric matrix. We thus obtain the following re¬ 
lation between parameter derivatives: 

[HuNi] p d = [N 2 ,H 2 ] P d- (4.43) 


Using (4.11), it can be rewritten as 



which shows that the parameter derivatives P\ p , p\d, Pip, and 
H 2 d are, in fact, not independent. The identity between pa¬ 
rameter derivatives similar to Eq. (4.43) was also reported in 
Appendix A of Ref. 38. 

By solving the eigenvalue problem (4.39) up to leading or¬ 
der for /rirf and fi 2 d, we obtain the following result: 

The dispersion relation and eigenvector corresponding to 
the Bogoliubov phonon are given by 

e = c p hk + 0(k 2 ), (4.45) 

( = w p h + ZphCphA: + 0(k 2 ) (4.46) 


D. Two-state approximation — quadratic dispersion 

In the above naive perturbation method, we cannot obtain 
the dispersion relations of ripplons in the finite-size system 
e ~ VZk 1 . In this subsection, we give a little better treat¬ 
ment to derive this. If the eigenenergy of the Bogoliubov 
equation is sufficiently small, only ripplon excitations exist. 
So, the eigenvector is well approximated by a linear combi¬ 
nation of two vectors, w T i p and z np . Using this fact, we solve 
the Bogoliubov equation non-perturbatively under the approx¬ 
imation such that the state space is spanned only by these two 
vectors. The result contains not only the previous complex- 
coefficient linear dispersion but also the yLkr behavior. How¬ 
ever, even in this treatment, we cannot obtain the dispersion 
relation and the eigenvector allowing to take the limit L —» oo. 
The final goal is given in the next subsection. 

Let us solve the Bogoliubov equation 

(H 0 + Mok 2 )^ = ^ (4-54) 

with the assumption that the eigenstate is given by the linear 
combination of the above two vectors: 

£ = awop +/hd- (4.55) 


with 



gnPi 

m i 



+ 


822P2 L _ A 
m 2 \ Lj 


w p h = VgiTwi + VI 22 W 2 , z p h = \z p . 


(4.47) 

(4.48) 


Here, “ph” means the phonon. Strictly speaking, the first order 
eigenvector z p h may include z,j, but we ignore it because it is 
not important in order for the first-order equation Ho£\ = ei£o 
to be satisfied up to 0(p l£ /). 

The dispersion relation and eigenvector corresponding to 
ripplons are given by 

e = Cn V k + Oik 1 ), (4.49) 

C, = Wrip + ZdpCripk + 0(k 2 ) (4.50) 


with 


c 2 _ (L 3 - d 2 )(pipu - PlPld) +of 2, 

Cnp mipi (L- d) + m. 2 p 2 (L + d) ,d 

Wn V = |l - m\Wi - |l + ^jm 2 w 2 , 

L 2 -d 2 

Zrip — l i Zd + k) (iZj) ■ 


(4.51) 

(4.52) 

(4.53) 


Since pipw - p 2 p 2 d is exponentially small and negative 
[Eq. (4.32)], c r ip is pure imaginary. Therefore, this dispersion 
relation represents the existence of unstable modes in a very 
narrow wavenumber region. 


Here, we use Zd instead of z np as a basis vector, since z r ; p oc Zd 
up to leading order with respect top,d’s [Eq. (4.53)]. Different 
from the previous subsection, the coefficients a and /j are now 
/.'-dependent. Taking the cr-inner product between Eq. (4.54) 
and w r ip, Zd, we obtain the 2x2 matrix equation 


-e 

1^2 (Wrip>^OWrip)<r 
V (}VripiZd)cr 


L % , i ,2 (ZdMgzj), N 
L 2 -d 2 (Wrip,Zd)a 


P 


= 0, 


(4.56) 


where we have used (z,-, Mowj)^ - 0, (w,-, MqWj)^ = d/j— for 

Lc^ 

i, j — 1 ? 2 3.nd HqZtvp — ^rip HqZ d — L 2 —d 2 • Let US 

introduce the notation 


T 0 := 


(■ Zd, MoZd)cr 
2 



(4.57) 


which represents the kinetic energy of the DW. By virtue of 
Eq. (4.12), the ^/-derivative takes up only the gradient energy 
of the DW, and it ignores the gradient energy near the bound¬ 
aries x — ±L. This means that the leading value of (zd, MoZd)cr 
does not depend on a choice of the BC for sufficiently large L, 
and hence it can be approximated by the kinetic energy of the 
DW in the infinite-size system: 



where we should consider t/q and t// 2 of the infinite-size system 
when we use Eq. (4.58). Then, solving Eq. (4.56) yields the 
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dispersion relation and the eigenvector 


e 2 = A 0 k 4 + c 2 p k 2 = A 0 k 2 (k 2 - k 2 c ). 

(4.59) 

L 2 - d 2 2 


£ ~ 1 & <-d-> 

„ . 2 (L 2 - d 2 )To 

(4.60) 

(4.61) 

m\p\(L - d) + ni 2 P 2 {L + d) 

k c := J c 2 /A 0 = 

\ V Y 2 T 0 

(4.62) 


respectively. Note that Aq = 0(L ) and k c is positive and of 
order 0 ( V7~e aL! 2 k) as a result of Eq. (4.32). The very narrow 
region 0 < k < k c gives the unstable modes. If the physical 
parameters of i/q and i //2 are symmetric, i.e., m\ = m 2 , d = 
0, pi = P2, p 2 d = ~Hid, it reduces to 


where /fow 2 is defined in Eq. (A18). We thus set 112 - V 2 = 0. 
As for F 1 and G 1 , if we assume (Fi,Gi) oc e ±lx , we obtain 
l = K\ + 0(Lr l ) and l - k + 0{L~ 2 ), where k\ is defined in 
Eq. (A12). The former solution is ignorable. The latter solu¬ 
tion can contribute and the corresponding approximate eigen¬ 
vector is given by 



/1 + 0(L~ 2 )\ 

i o(L-y 2 ) j 


n ±kx 


(4.68) 


Thus, we can set F 1 = e ±kx and Gi = 0, implying that the den¬ 
sity fluctuation is ignorable, as stated in the procedure (B) of 
Sec. IIC. Moreover, following the procedure (B), we impose 
the Neumann BC at x — -L. Then, we have 

Mie~ lSl = -vie 1 ® 1 = coshk(x + L ), U 2 - v 2 - 0 (4.69) 


LT °-k 2 (lc - k 2 ), k c = 


mpi 


-POP Id 


To 


(4.63) 


If the narrow complex region is ignored, it gives e = J -^-k 2 . 


mipi 


as Ref. 32. (Note that the mass is taken as 2mi = 1 in Ref. 32.) 


E. /.-dependent two-state approximation — interpolating 
formula 


The approximations used so far could not produce disper¬ 
sion relations and eigenvectors which allow to take the limit 
L —> 00 . To accomplish this, let us construct a modified quasi¬ 
particle wavefunction including the asymptotic behavior far 
from the DW, corresponding to the general procedure (B) in 
Sec. IIC. Let us consider a uniform region —L + £ < x < d — ij 
so that the approximate expression i/q = \fp\t l6 ' = const, 
and tj /2 — 0 can be well applied. Here £ is a typical healing 
length of the condensates. We further introduce the notation 
F\ = iqe -1 ® 1 - vie 1 ® 1 , Gi = jqe -1 ® 1 + vie 1 ® 1 . Then, in this 
uniform region, the Bogoliubov equation can be written ap¬ 
proximately as 


-d 2 x + k 2 

2 m 1 


-F 1 = eGi, 


~d 2 x + k 2 

2m 1 


+ 4gnpi Gi = eF 1 , 


-d\ + k 2 

2l1l2 

-d\ + k 2 
2/722 


+ 2(gi2Pl - g22P2)J "2 = eu 2 , 
+ 2(gl2Pl - g 22 P 2 ) 1'2 = ev 2 . 


(4.64) 

(4.65) 

(4.66) 

(4.67) 


Let us find a solution under the approximation such that we 
ignore functions whose decay rates are comparable with the 
healing lengths of condensates. (See Appendix A for ex¬ 
pressions of the healing lengths.) We are interested in the 
wavenumber of order k ~ 0(Lr l ). Correspondingly we as¬ 
sume e ~ VZk 2 ~ ()(L !,/2 ). In this approximation, U 2 and 
V 2 are ignorable, because if we consider the solution 112 , V 2 °c 
e ±lx , we obtain 1 — (kq W2 + k 2 ± 2 m 2 e )^ 2 = kdw 2 + G(Zr 3 ^ 2 ), 


for the region -L + £:< x <d - £. By the same argument, 
in the right-side uniform region d + £<x < L - £, assuming 
i/q = 0 and \jj 2 = -^02 e 1 ® 2 , we obtain 

Mi = vi = 0, M 2 e^‘® 2 = -V 2 e‘® 2 = coshk(x- L). (4.70) 

We thus obtain 


'u 1 ' 


' a 6 {d - x)e lSl cosh k(x + L) ' 

U 2 


b 0 (x - (i)e 1 ® 2 cosh k(x - L ) 

Vl 


-a 6 {d - x)e -1 ® 1 cosh k(x + L) 

,V 2 , 


-b 0 (x - d)er iei cosh k(x - L), 


where the coefficients a, b are fixed below. This conclusion 
is more quickly obtained if we assume that e is small hence 
ignorable. 

Next, by following the procedure (C), we modify the solu¬ 
tion (4.71) to include the zero-mode solution Wri p [Eq. (4.52)]. 
Henceforth we write such modified solution as wo P (k). The 
modified solution must satisfy Wri p (0) = Wnp. From the 
expression (4.71), we can conceive the replacement 6 (d - 
x)e lBl —2 <Ai / y[p\, 0 (x - fif)e‘® 2 —> ^ 2 / ^Jp 2 to include Wnp. 
Then, we obtain 


' a’tpi coshk(x + L) ' 
b'lfj 2 cosh k(x - L ) 
-a'l/s* coshAfx + L) 
-b'xf/ j coshAfx - L), 


(4.72) 


Here a' = a / ^[p\ and b' = b/ ^[pF. The ratio of the coeffi¬ 
cients a', b' is fixed by imposing the condition that w r i p (k) has 
the same behavior with Wri P near the DW, that is, 


H'rfp(k) - wnp for x ^ d. (4.73) 


Then, we have 


tVrip(^) — 


(4.74) 
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It is worth noting that this solution can be used for both 
Dirichlet and Neumann BCs. Now we solve the Bogoliubov 
equation by the modified ansatz 


Here Ao is introduced in Eq. (4.61) and its size-dependence 
is Ao = O(L). Correspondingly, the dispersion relation (4.77) 
reduces to 


C = awn V {k) +fizd- (4.75) 


If we set A = 0, i.e., Wri p (A) = Wri P (0) = Wri p , the ansatz 
reduces to that in the previous subsection [Eq. (4.55)]. By 
taking the cr-inner product between the Bogoliubov equation 
(Hq + Mok 2 )£ = and w r i P (0) and Zd, we obtain 


ffrig (~rf.M’rip(0))ir 7.2 (ZdMgZd)cr 

L 2 -tP (Zi.WnpCtJV (H.WhpWV 

i 2 (Wrip(0),MoWrip (&))<r 

i. (WripfOi.ZiV € 

(4.76) 



where we have used the easily-checked relations 
(Wri p (0), WripWV = (.Zd, M 0 Wn P (k))a- = 0. The disper¬ 
sion becomes 


e 2 = A(k)lc(k 2 - A 2 ), 

2r 0 (w r i P (0),M 0 w r i P (A)) o - 
k) :=---, 

(Zdi VV, i p (0) ) fr (2^, VV’i j p (k)) (r 


(4.77) 

(4.78) 


e 


2 


A 0 A 2 (A 2 - A 2 ) 
2T 0 k 3 

m\p\ + 1U2P2 


(kL <sc 1) 
(L —> oo). 


(4.85) 


We thus have found that the dispersion relation (4.77) includes 
both e ~ y/Lk 2 and e ~ A 3,/2 . Furthermore, the formula (4.77) 
with Eq. (4.83) is valid even for the intermediate region inter¬ 
polating these two limiting cases. 

If the DW is located at the center (d = 0), the expression 
for the dispersion relation becomes a little simpler: 


27q tanhAL , , , 

- ~ - -k 2 (k 2 -k 2 y 

m i p i + m 2 p 2 k 


(4.86) 


It includes all three behaviors shown in Table I: 


2T 0 


I mipi + m 2 p2 


x 


i VZk c k, 
VIA 2 , 

k 312 . 


0 < k < 0(e- aLl t), 
0(e- aL/ f) <k< 0(L~ l ), 
0(L- X ) < k < 0(r l l 

(4.87) 


where k c is defined in Eq. (4.62). Let us evaluate the leading 
order of A-dependent cr-inner products appearing in A (A). In 
fact, the following rough expression is sufficient for this pur¬ 
pose: 


|^i| 2 = Pl 6(x + L)6(d - x), (4.79) 

I 1 A 2 I 2 = PiS(x - d)6(L - x). (4.80) 

We emphasize that these expressions should not be used to 
evaluate other cr-inner products such as 2Tq = (z.d, M( } Zd), T - 
By using Eqs. (4.79) and (4.80), we obtain after some calcu¬ 
lations: 


(.Zdi W'ii p (A)) 0 - — (z.d , W’npCO )) (r 


= Wipi(l - f) + OT 2 p 2 (l + f), (4.81) 

(wrip(0), M 0 Wrip(A)) o . 


d 2 tanhA(L + c/) d 2 tanh k(L - d) 

= wipi(l - f) ---+ m 2 p 2 (l + f) ---- 

(4.82) 


We thus obtain 
A(k) = 


2T 0 


k[m\p\(L - d) + m 2 p 2 (L + d)] 2 
|/nipi (L - c/) 2 tanhA(L + d) 

+m 2 P 2 (L + d) 2 tanh A(L - rf)]. (4.83) 

This A(A) has the following two important limiting cases: 




A(A) = 


2T 0 


y k(ni\p\ +m 2 p 2 ) 


(AL«c 1) 
(L —> 00 ). 


(4.84) 


The eigenvector is given by 

A 2 

£ = ewrip(A) + 


k-A(k)[m l p l (l - f) + m 2 p2(l + ()] 
2To 


-Zd- 

(4.88) 


If we set d = 0 and take the limit L —> 00 , we obtain 






1 


d 

dx 

^2 

r, 

/ 2T ° 7 , 1/2 

y m\p\ + m 2 p2 

-m 2 i)/2^ kx 

-mn//\e kx 


w 2 ) 


m 2 i// * 2 e~ kx 


(4.89) 


where Eq. (4.33) is used. It describes the quasiparticle wave- 
function of ripplons in the infinite system. 


V. SUMMARY 


In this paper, we have presented the analytical formulas in¬ 
terpolating the integer dispersion in finite-size systems and 
non-integer dispersion in infinite-size systems for the Kelvin 
modes along a quantized vortex and the ripplons on a domain 
wall in superfluids, together with quasiparticle wavefunctions, 
and have found a complete agreement between our formulas 
and numerical simulations. The derivations of these formulas 
are supported in a fully analytical way using the techniques 
constructed in Ref. 32. 

Finally we give a remark on the criteria for emergence of 
non-integer dispersion relations. Inferromagnets,NGMs such 
as a ripplon on a domain wall 39 and Kelvon on a skyrmion 
line 40,41 have quadratic dispersion relations even for large sys¬ 
tem sizes. This is because the zero modes in these systems are 
normalizable. On the other hand, in the cases studied in this 
paper, the zero modes are non-normalizable 32 . 
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Appendix A: Healing lengths of two-component BECs 


In this appendix we discuss a few fundamental facts on 
the two-component BEC model such as conservation laws 
and healing lengths of the DWs. Let us consider an infinite 
one-dimensional system. The time-dependent GP equation is 
given by 

= (-pi - ^d 2 + 2gnW\V + ^gnWiV^x, (Al) 

id t l //2 - (~p 2 - + 2*211^1 r + 2 .g 22 \^ 2 \~)^ 2 - (A2) 

Here we write down the conservation laws. The number con¬ 
servation laws are 


l2 . ^ (iWjlfi-Wix) 

a 1 * 1 + at(-St- 


= 0, 


i = 1,2. 


(A3) 


The momentum conservation law is given by 


d_ 

dt 


Z 

>1,2 

d 


i oioi'i - r^u) 


dx 


z 

i=l,2 


iWa-WO . Mix] 2 . , , ,2 

-~-+ -z — +p#;' 

Z Zrrii 


- 2 gyWV/ 

i,7= 1,2 


= o. 


(A4) 


We omit the energy conservation law because it does not give 
a new integration constant for a time-independent solution. 
From these conservation laws, for the stationary solution i/q, = 
ip 2 t = 0, we have the following integration constants: 


. i (Tuti-Wix) . , „ 
Ji= ---, i= 1,2, 


jmom ^ ^ 
i=l,2 


2/72/ 
\^ix\ 2 


(A5) 


2/77/ 


■/til^.i 2 ]- gyWV/- (A6) 

' i,j= 1,2 


If i [/\, i [/2 are real, j\ = j 2 = 0, and hence y mom is the only 
non-trivial constant. 

Let us consider the DW solution having the following 
asymptotic form: 

|° (x^+oo) (x^ + oo) 

\ #r (*-»-«>), V2 }o (x > —oo). 

(A7) 

In order for this asymptotic form to become the solution of the 
GP equation, the values of the chemical potentials should be 
fixed as 


Hi = 2g u pi, i = 1,2. (A8) 

Furthermore, from the x-independence of the momentum cur¬ 
rent density (A6), we obtain the relation 

jmom — gllPl — g 22P 2 , (A9) 

which is the same with Eq. (4.9). Thus, p\ and p 2 cannot be 
chosen independently. We also note that the meaning of the 
parameter p is, in fact, the square root of the momentum cur¬ 
rent: j mom = p 2 . 

Let us introduce four kinds of healing lengths. We first 
consider the situation such that only i/q exists. In this case 
Eq. (A6) reduces to 

*Aia = 2mign(^ - pi) 2 , (A10) 

and a solution is given by the well-known dark soliton solu¬ 
tion: 

*At = \4oTtanh^, (All) 

ki := 2y/2g n m 1 p l . (A12) 

This k\ describes the inverse of the healing length for the one- 
component system. In the same way, we obtain that for xjtp. 

K 2 := 2 xj2g 2 2m 2 p2- (A13) 

Next let us consider the decay rate of i/q on the right side of 
the DW, where t/q is dominant. Assuming t//| is small and 
\j /2 — \[p 2 , the GP equation can be approximated as 

^4 + 2(g l2 p 2 - gnPiMi = 0, (A14) 

where the nonlinear term is ignored with assuming small i/q. 

Then, 


i/q cx e _ ' fDwlJC , 

/cdwi := 2 ^/77ip2(gi2 - yfgug 22 )- 


(A15) 

(A16) 


This /fowl represents the decay rate. Here, we have used 
Eq. (A9) to obtain g n p 2 - gupi = P 2 (gn ~ y/gngn)- By 


the same calculation, on the left side of the DW, we can show 
ij /2 oc e* DW2 *, (A 17) 

*"DW 2 := 2. Jm 2 pi(gi 2 - yjgng 22 )- (A18) 
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Summarizing, we have obtained four inverse healing lengths, 
Ki, K 2 , k DW i, and /q 2W2 - Thus, the term “typical healing length 
used in Secs. II and IV precisely means the largest one 
among these four lengths, i.e., 

f — max(A' 1 ,/e 2 i^dw 1 ’^dw 2 )- (A19) 


i f/{ occupies is given by 


Li = aK(mi)yj^-, 

(B7) 

(1 (the Neumann BC) 

12 (the Dirichlet BC). 

(B8) 


Appendix B: Evaluation of k c for the case of gn = 00 

In this appendix, we focus on the system with g 12 = 
+ 00 , in which two condensates are completely decou¬ 

pled and hence the GP equation reduces that of a single¬ 
component BEC. We want to find the leading L-dependence of 
k c [Eq. (4.62)], the maximum wavenumber such that the dis¬ 
persion relation of ripplons becomes complex-valued, in other 
words, the maximum wavenumber of unstable modes. For 
simplicity, we only concentrate on the case where the physi¬ 
cal parameters of two BECs are symmetric, i. e., gn = g 2 2 = 
1, 2m\ - 2/722 = 1, N\ = N 2 ■ In this case, pu — -p 2 d holds 
by symmetry. 

Both i/p and iA 2 satisfy the single-component GP equation 
—if/" -/rtA + 2 if? = 0, (Bl) 


and the general solution is given by 


ils(x-p,m)= 

rnj , (B2) 

E(m) 

Q(m) = 1 “ K{mY 

(B3) 

(1 + m)p 

^ Q(m) 

(B4) 


Needless to say, L\ and Li are not independent and satisfy 
L\ + L 2 — 2 L. 

Since we want to solve the energy minimization problem 
with respect to L\ under the condition that N\ , Vi, L are fixed, 
we change the independent variables from m\,pi,m 2 ,p 2 to 
Ni,Li,N 2 ,L 2 . Their relations are given by 



^ = KUmfQOn,). (BIO) 

Thus, in order to move on to the description by V, and L,, 
we need an inverse function of K(m) 2 Q{m). Though the exact 
inverse function cannot be written down in a closed form, if 
m - 1 (i.e., if sn is almost tanh), we obtain the following 
asymptotic expansion: 

x = K(m) 2 Q(m) (B11) 

<-» m = 1 - 16e - - v + 128e~ 2y + • • • , y := 1 + Vl + 4.r. 

(B12) 

The expansion (B12) can be obtained by using the formulas 

K( 1 - 165) = - ^ log 5 - 25(2 + log 5) + 0(6 2 log 5), (B13) 
£(1-165)= 1 - 45(1 + log 5) + 0(6 2 log 5) (B14) 


where K and E are the complete elliptic integral of the first and 
second kind, respectively. Here and hereafter, we use Math- 
ematica’s notations for the elliptic integrals/functions unless 
otherwise noted. The solution (B2) is characterized by two 
parameters m and p. The former is an elliptic parameter and 
satisfy 0 < m < 1. The latter has the physical meaning of the 
averaged particle number density: 


f 


'K(m) y Q(m)/p 


K(m) ylQ(m)/p 
The energy per particle can be calculated as 

rK(m) 2jQ(m)/p 

Jo 


dj#|“ = p. 


(B5) 


E 

N 


.'|2 . 


fo 


dx(|l/r' 


) [in + (1 + m)Q(m)]p 


3 Q(m ) 2 


(B6) 


Henceforth, we write the physical parameters of ip\ s (i =1,2) 
as mi, pi, pi, Nj, Ei, and so on. 

If we use the Dirichlet BC (<//,■ = 0 at the boundary), the 
profiles of i//j’s are given by the sn function with one-half of 
a period. If we use the Neumann BC (t^' = 0 at the bound¬ 
ary), the profiles of i/r, ’s are given by the sn function with one- 
quarter of a period. Therefore, the length L, of the region that 


and solving the equation x = K 2 Q = K 2 - KE w.r.t 5 iter¬ 
atively. When Eq. (B12) is applicable, K(m) and Q{m) are 
given by 

K(m) = y Q + 2e- y - 16e^ 2y + ■■■), (B15) 

c<m) = iSF = (‘"9( 1 “ 8e ‘’ +112e ‘ 2 ' + "T 

(B16) 


By using them, the chemical potential and the energy for 1 In 
are written as a function of (L,, Nj): 


11 

(l + r) (l _ 48e~ 2y + •••), 

(B17) 

N 2 r 

£ ' = 77l 

(>+5)-?( 4+ 7k 2 T 

(B18) 

y : — 1 + 


(B19) 

Here, the terms of order 0(y a e 2hy ) with a > 2 or b > 2 are 


ignored. 

Now, let us write 


Li = L + 5L, L 2 = L - 5L, (B20) 

Vi = Lp 0 + 6 N, N 2 = Lp 0 - 6 N. (B21) 
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where po = N '^ 2 is the average of the total particle number It is obviously negative: p\,j < 0. By ignoring the 0(1) nu- 

density. Let us minimize merical factor, the main L-dependence can be given by 


Ltotai - Ei + E 2 (B22) 


with respect to 8 L under the constraint that L,po, and 8 N are 
fixed. If 8 N = 0 <-» N\ = N 2 , we immediately obtain a trivial 
solution 5L - 0. Let us find 5L for the non-zero imbalance 
8 N 4 0. After a little tedious calculation, we obtain 


((Llotal 

85L 


<-> 


SL^SN 




lQ24L 2 ppe~ 2 ' 


4Ly/P0 


3a 2 


c 


9 a ) 


+ 


0(e“ 


8 L^JPO 




(B23) 


where 0{L 2 ) terms are ignored in each parenthesis. By using 
this 8 L, up to the same approximation, pi can be written as 

„ 4L^q 

* = + ife) - SN 2mp %- ° (1 - j^=). (B24) 

In the present calculation, recalling that we have set gn = 
g 22 = 1, the parameters p and cl introduced in Subsec. IV A 
are 


5N 

p =p 0 , d = 

2po 


(B25) 


Thus, the ^-derivative of p\ up to leading order is given by 

(B26) 


P\d 


dpi 

dd 


4096 pzL 

-e 


iL 4P0 


3 a 2 


J-Le (a = 1; Neumann BC) 

| -Le- 2L ^ (a = 2; Dirichlet BC). 


(B27) 


Since k c oc y-p w [Eq. (4.63)], we also obtain 

j VLe~ 2L ~^" (a = 1; Neumann BC) /-R 9 o\ 
e ~ | VLe- L ^ (a = 2; Dirichlet BC). 

We thus have proved the behavior in Eq. (4.32). 

Though this result is rigorously applicable only for the spe¬ 
cial case g 12 = oo, the numerical results suggest that the above 
behavior is also true for finite g\ 2 if we modify the exponential 

21-y^Q _ 2Ly/pQ 

factor as e « —> e v « , where v ~ 1 is a numerical fitting 

parameter. See Fig. 12. Thus, we can say that k c is always 
exponentially small. 

The above result suggests that the Neumann BC can sup¬ 
press unstable modes more strongly than the Dirichlet BC. 
For example, if we set L = 12 and po = 1, then k c ~ 10 7 for 
the Dirichlet BC and k c ~ 1 (L 10 for the Neumann BC. This 
means that the typical eigenenergies of the complex-valued 
regions are given by |e[ ~ Oik 2 ) ~ 10 1(1 for the Dirichlet 
BC and |e| ~ 0(k 2 ) ~ 10 2l) for the Neumann BC. While the 
former might be numerically seen, the latter is impossible to 
detect in the usual precision. Therefore, the Neumann BC is 
a powerful tool if one is interested in the infinite-size physics 
and wants to ignore finite-size effects, though sometimes this 
BC is not physically realistic. This observation is consistent 
with the previous numerical study performed in the Neumann 
BC in Ref. 18, where no unstable mode was found numeri¬ 
cally for large L. 
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